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ABSTRACT 



We investigate the evolution of a massive black hole pair under the action of dynamical 
friction by a uniform background of light stars with isotropic velocity distribution. 
In our scenario, the primary black hole Mi sits, at rest, in the center of the spherical 
star distribution and the secondary less massive companion M2 moves along bound 
orbits determined by the background gravitational field. M2 loses energy and angular 
momentum by dynamical friction, on a time scale longer than the orbital period. 
The uniform star core has total mass Mc and radius Tc, and the following inequality 
Me > Ml > M2 holds. 

In this paper, we investigate mostly analytically the secular evolution of the orbital 
parameters, and find that angular momentum (J) and energy {E) are lost so as to 
cause the increase of the eccentricity e with time, during the orbital decay of M2. 

In the region of the core where the motion of M2 is determined by the mean 
field generated of the uniform stellar distribution, E and J are lost exponeniially 
on a time scale ~ tdf determined by the properties of the ambient stars. The rise 
of e establishes instead on a longer time ~ tdf {Tc/rxY increasing as the apocenter 
distance ta decreases. 

With the progressive decay of the orbit, M2 enters the region r < ~ (Afi/Afe)^^'^^c, 
where the gravitational field of the primary black hole dominates, but the star back- 
ground maintains still uniform (to first approximation). Inside r^, a key parameter of 
the calculation is the ratio between the black hole velocity v and the stellar dispersion 
velocity a. (i) If v < cr, energy and angular momentum are lost exponentially on 
a time scale tdf- The growth of e occurs instead on a time r*^ longer than tdf by 
a factor ~ [a/vY'. Therefore, e rises weakly during orbital decay, r*^ is also found 
to be a function of e and increases as e ^ 1. (ii) In the opposite limit, i.e., when 
f > cr, the evolution of E and J is close to a power-law and establishes on a time scale 
~ [v/aY ty,f- The eccentricity grows on a time comparable to this scale. Along an 
evolutionary path, e increases significantly: This rise leads the pericenter distance to 
diminish exponentially, in this limit, t'^ is a function of e and decreases as e — ^ 1. 
This limit {v > a) is attained close to the cusp radius rcusp ~ {Mi/Mc) r,.. i.e., the 
distance below which the stellar distribution is affected by the gravitational field of 
Ml. Below Tcusp our description is invalid, and we terminate our analysis. 

Energy losses by gravitational wave emission become comparable to those by dy- 
namical friction at a critical distance that depends on the ratio Mi/Mc- Consistency 
with the model assumptions implies Mi -C Mc. The braking index n is calculated 
in this transition region: A measurable deviation from the value of 11/3 correspond- 
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ing to pure gravitational wave losses provides ideally an indirect way for probing the 
ambient medium. 
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1. INTRODUCTION 



Galaxies in their cores may contain massive binary black holes. This hypothesis 
follows mainly from two considerations, (i) Most galaxies harbor a central black 
hole, relic of an earlier active phase (Dressier 1989; Bland-Hawthorn, Wilson & TuUy 
1991). This assertion is plausible if the AGN phenomenon is the result of short- 
lived phases of activity that have involved most of the galaxies in the past (Rees 
1990). (ii) Galaxies interact and merge. Evidence comes from the appearance of 
unusual morphological features interpreted as clear signature of a tidal interaction 
between companion galaxies or of a recent merger (see Barnes & Hernquits 1992 for 
a complete review). If galaxies participate a hierarchical process of clustering, large 
galaxies cannibalize the small ones and in these systems a black hole binary with 
unequal masses may form (Begelman, Blandford & Rees 1980, BBR hereafter; Roos 
1981; Gaskell 1985). Recently, a claim was put forward for the observational evidence 
of a massive binary system in the nucleus of 1928+738 (Roos, Kaastra & Hummel 
1993). 

It is believed that the nuclei of two interacting galaxies sink, during a merger, to- 
ward the center of the common potential well by dynamical friction: a gravitationally 
bound pair of black holes may thus form, embedded in core stars (BBR; Valtaoja, Val- 
tonen & Byrd 1989; Governato, Colpi & Maraschi 1993). The frictional deceleration 
on the pair itself progressively tightens the system and the two black holes eventually 
bind. When their distance becomes comparable to the mean star separation, single 
interactions with stars on a time less than the orbital period come into play and en- 
ergy is transferred to the light particle. This process however may lead to a depletion 
of background stars that are driven into loss cone orbits via a sling-shot mechanism 
(BBR; Roos 1981, 1988). The decay of the binary may be relented or even inhib- 
ited at separations where the energy loss by emission of gravitational waves would 
be weak to drive the binary to coalescence within the Hubble time. Mechanisms for 
replenishing the loss cone are therefore invoked for triggering further evolution (BBR; 
Ross 1988; Polnarev & Rees 1994). 

Recent numerical iV-body simulations on the dynamical evolution of black hole 
binaries in the violently relaxed core of a merged galaxy (Fukushigc, Ebisuzaky & 
Makino 1992a; Mikkola & Vahonen 1992; Makino et al. 1993), have indicated that the 
collective effect of dynamical friction leads to a "runaway" increase of the eccentricity 
of the binary, and to a steep drop of the pericenter distance. Losses by gravitational 
wave emission would thus intervene early in the life-time of the binary that evolves 
rapidly (within 10^ yrs) toward final coalescence. Proposing an approximate extension 
of the Chandrasekhar formula for dynamical friction in a non-uniform anisotropic 
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medium, Polnarev & Rees (1994; PR hereafter) studied numerically and analytically 
the evolution of a black hole pair in a typical galaxy. They find that in the core 
regions where sufficiently large density gradients develop the eccentricity decreases 
and velocity anisotropics tend to stabilize this tendency. 

The study of the black hole binary evolution is important since a correct deter- 
mination of the life-time and occurrence rate of coalescence (Fukushige, Ebisuzaky 
& Makino 1992b) is crucial for assessing the relevance of black hole binaries as can- 
didate sources for the detection of gravitational waves (Thorne & Braginsky 1976; 
Thorne 1992). These sources could indeed be detected by Doppler experiments by 
spacecraft tracking (Estabrook & Wahlquist 1975) and long interferometer in space 
as LISA (Danzmann et al. 1993) and SAGGITTARIUS (Hellings et al. 1993). 

Moving from these considerations, we here re-explore the orbital evolution of a 
black hole pair, using mainly an analytical technique. In our model, a massive black 
hole Ml sits, at rest, in the center of a uniform galactic core with total mass Mc, 
and radius Tc (Binney & Tremaine 1987). A comparatively low mass black hole M2 
moves along bound orbits spiraling inward under the action of dynamical friction by 
the field stars. We consider the relevant case when M2 < Mi < Mc and introduce 
the dimensionless parameter 77 = Mi/Mc < 1. 

Along the course of dynamical evolution, the frictional deceleration on M2 can 
be considered as small perturbation relative to the "background" gravitational field, 
and therefore causes only secular changes in the orbital motion. In our scenario, 
M2 enters the core region of radius Tc along bound orbits as it experiences the mean 
gravitational field generated by the uniform stellar distribution (see PR for detail). 
M2 eventually bind to the massive central black hole, at a distance ~ r/^/^ r^. 
Below this separation, M2 is still surrounded by the homogeneous sea of field stars. 
The influence of the central black hole on the stellar distribution becomes important 
at a smaller separation rcusp ~ i^cV- Below rcusp steep stellar density gradients and 
anisotropics in the velocity field appear, and our treatment ceases to be valid. 

In this paper we determine the evolution of the orbital parameters of M2 in the 
two relevant regions: (a) the "outer" zone, in the interval rB < r < rc, where the 
mean gravitational field is generated by the core stars and (ii) the "inner" zone 
'"cusp < r < rB where the gravitational field of the massive central black hole affects 
the motion of M2. 

The plan of the paper is as follows. In §2 the expression for dynamical friction 
is reviewed and key parameters of the model are introduced. In §3 we outline the 
basic equations for the energy and angular momentum losses. We solve these equa- 
tions using Gauss' perturbation theory, both in the "outer" and "inner" core regions. 
Characteristic time scales are derived as a function of the eccentricity. A criterion 
for the growth of e is also presented and applied to our analysis. In §4 we describe 
possible evolutionary paths the binary experiences, within the inner region of the 
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galactic core. In §5 we review formulae on the emission of gravitational waves: the 
coordinated effect of dynamical friction and emission of gravity waves is then dis- 
cussed. The deviation, due to dynamical friction, of the braking index n from the 
canonical value of 11/3 is calculated. §6 contains our conclusions. 



2. THE MODEL 

The frictional deceleration experienced by the secondary BH in a homogeneous back- 
ground of stars with isotropic gaussian velocity field is 

Vdy = -xa{x)v , (2.1) 

where ■& is a unit vector in the direction of the velocity of the BH relative to the 
background. 

In equation ( p.l| ), the BH velocity is expressed in unit of the one-dimensional 
dispersion velocity (-\/2 a) of the stars 



x = ^ (2.2) 



and the function g{x) reads 



erf(x) -^xe ^'^ 

71 



9i^) = ^ 

(Chandrasekhar 1943; Binney & Tremaine 1987). 

The coefficient x contains key parameters of the background stars 

^ , G^mriQ (M2 + m) 
X = 27rlnA 



(2.3) 



_g / Hi ii \ ; M2 \ f m \ 1 300kms~^\^ _2 



where m (-C M2) is the mass of a typical individual star, (~ Mc/mr^) the number 
density of stars and InA = ln{bmax/bmin) the so called Coulomb logarithm (Binney & 
Tremaine 1987). 

The dissipative force (eq. [2.1]) parallel to the velocity vector v clearly acts so as 
to decelerate the secondary BH causing the progressive decay of the binary separation. 
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Equation (|2.1| ) takes a simple asymptotic expression in two relevant limits. In the 
first limit, i.e., for t> <^ a (i.e., x <^ 1), 

4 3 

Vj)F = -i^xix- -x-^ + o{x-^))v. (2.5) 
Sy'n 5 

This regime may characterize the early phase of the BH evolution when orbiting 
regions where field stars move with dispersion velocity a higher than v (BBR). 

In the second limit, i.e, for v ^ a (i.e., for x ^ 1), the deceleration takes the form 

■"DF^-X^^^- (2.6) 

The transition between the two asymptotic regimes occurs when a; ~ 1, corresponding 
to a value of the mean separation between the BHs ri ~ rjVc- 

Simulations have shown that equation (2.1) often provides an accurate descrip- 
tion of the deceleration experienced by a single object. In our model for the binary 
evolution, the drag acts only on the secondary black hole since the primary is at rest 
relative to the medium. If the primary were to move a correction to equation (2.1) 
would be of need (Bekenstein & Zamir 1990; see however Gould 1993). 

A source of uncertainty in equation (2.1) is the evaluation of the Coulomb loga- 
rithm, and in particular the estimate of the maximum impact parameter bmax whose 
value is predicted from numerical simulations to be in the range iV-^/3^c (Smith 
1992) and Vc (Farouki & Salpeter 1982): We have chosen as reference value InA = 10 
(obtained considering bmax = ^c] Binney & Tremaine 1987). The deviation related to 
a different choice of bmax (i-e., a fractional error of less that 30%) falls anyway within 
the intrinsic spread of InA related to the typical parameters of galaxy models. 

Equation (2.1) is considered valid down to the limiting radius r^usp — V'^c- At 
this distance the influence of the primary BH on the stellar distribution is important 
(Frank & Rees 1976; Binney & Tremaine 1987) and different mathematical tools are 
of need to follow the subsequent evolution (PR; see Maoz 1993 for a new approach to 
DF). 

2.2 The Secondary Black Hole Orbital Motion 

(a) The " outer" region: In absence of dissipative forces the secondary BH un- 
dergoes a periodic motion. In the region tb < r < rc, this motion is determined by 
the gravitational potential U{r) generated by the uniform and isotropic distribution 
of stars. According to the Poisson equation 

t/(r) = lfiV, (2.7) 
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where Q is related to the density of stars p — mno through the following relation 



n'^^Gp. (2.8) 

The secondary BH undergoes a periodic motion along elliptic orbits: In cartesian 
coordinates, with origin in the center of the star distribution, the motion is described 

by 

X = ta cos(Qi) , (2.9) 

?/ = rp sin(Qt) , (2.10) 

where ta and rp are the apocenter and pericenter distances, respectively. 

We can express the two constants of motion, i.e., the total energy E and angular 
momentum J in the following form 

= AST- ("2' 

and introduce the eccentricity 

e^'-^^. (2.13) 
TA + rp 

For the analysis of the motion of the secondary BH driven by DF it is of importance to 
estimate the magnitude of the deceleration -udf and compare it with the acceleration 
generated by the gravitational potential = fl'^r. In the limit a; <^ 1, we have 

^ = -^GM2\nAJG~pa-^ 

= 5 X 10-3 (^] ( C-^f' /300kms-^\- 

The weakness of the frictional drag enable us to treat DF as a perturbing force. 

For the harmonic potential, we consider only the case a; -C 1. From the virial 
theorem and from equation (2.8) this condition is equivalent to 

— < 1 (2.15) 

which is obviously satisfied. The opposite limit {x ^ 1), equivalent to r-p/rc > 1, 
corresponds to a motion confined in regions that are external to the galaxy core. This 
limit is therefore inconsistent. 
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(b) The "inner" region: The secondary BH becomes eventually bound to the 
primary BH. This occurs at a separation 



\10^Mq ) UOOpc^ 



TB^r]l re ~ 37 — — — — pc (2.16) 



for typical parameters of a galactic core (BBR). At radii r < r-Q, the potential of 
the central BH dominates over the mean potential generated by the homogeneous see 
of stars. At rs the BH binary forms and the secondary BH moves along Keplerian 
orbits. In the region rcusp < r < r^, the binary is surrounded by ambient stars having 
a uniform density and dispersion velocity cr larger than v. This is the reason why 
only below the cusp radius 

the stellar distribution is affected by the central BH. 
The total newtonian energy of the binary 

E^y-Et (,18) 
la 

and angular momentum J of magnitude 

= GMii^a{l - e^) (2.19) 

are constants of motion (in eq. [2.18] and [2.19] M = Mi + M2 ~ Mi and /i ^ is 
the reduced mass). Given E and J the semimajor axis a and the eccentricity e are 
uniquely determined. 

Analogously to the previous case, we estimate the magnitude of the frictional drag 
■Udf relative to the gravitational acceleration, v^^ — GM/a^. For x <^ 1 we have 



37rV2 a {GMyi'' 



^''^ ^300kms-i\^ /lO^M©^^ 



and for x ^ 1 



VK " (CM) 
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Even in this case DF is weak compared to gravity. Therefore, DF can be considered 
as a perturbing force that controls only the evolution of the BH on a time scale 
longer than a orbital period. In both potentials, the motion of the secondary BH 
is almost periodic. On the secular scale fixed by DF, the mean motion varies, the 
system evolving toward states that are progressively more bound. 



3. THE LONG TERM EVOLUTION OF THE SECONDARY BH 

3.1 Long Term Variation of E and J 

In this section we derive the long term evolution of E and J. In presence of DF, energy 
is not conserved and it is dissipated at a rate 

E — fi vdf ■ V 

= -f^X9ix)v. (3.1) 

Similarly, the equation for the angular momentum loss is 

j = iJ,{r X vdf) 

= -X9{x)-. (3.2) 

V 

E and J represent the instantaneous values of the time derivatives for E and J. 
Since DF is a weak perturbing force, it is meaningful to calculate (E) and (J), i.e., 
the average on the orbital period. This mean is defined for any orbital parameter S 
(such as E, J, rp, ta, a and e) as 

{S) = ]^rSdt, (3.3) 
Jo 

where P is the orbital period. Below, the results are presented separately in the two 
asymptotic limits of equation (2.1) and for the two expressions of the gravitational 
potential acting on the secondary BH: the calculations are fully analytical. 



(a) The Harmonic Potential 



When the gravitational potential of the star distribution is dominant (r^ < r < 
Tc) , the losses can be evaluated for x <^ 1 retaining only the first term in the expansion 
of g{x). This yields: 
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It is relevant to notice that in this hmit the decay is exponential as DF is proportional 
to (f ) in this limit (see eqs. [2.5], [3.1] and [3.2]). These losses are also to lowest order 
independent of the eccentricity e. 

As indicated by equations (3.4) and (3.5), the magnitude of the decay time is 
primarily determined by the values of x a-nd a. We are therefore led to introduce a 
characteristic time scale for DF 



tdf 

\ 

/ ^ \3/in6/\^\/owin3 ^„-3\ 

yrs . (3.6) 

In terms of tdf the decay time scales for E and J are equal and read: 

r| = ri = :^rDP . (3.7) 



(b) The Keplerian Potential 

Below tb the gravitational potential of the massive primary BH dominates and 
the motion of the secondary is nearly Keplerian. 

(i) The limit x < 1 

We consider the case of a; ^ 1 and, for simplicity, retain again the first term in 
the expansion of g{x). In this case the losses of energy and angular momentum read 

(the details are presented in Appendix A and Appendix B. In the calculations x is 
a function of the orbital phase and the asymptotic limit is meant to be valid over 
the complete orbit). It is important to notice that even in this case the losses are 
exponential and independent of the eccentricity. According to the previous equations 



11 



the corresponding time scales are 



.E 



-T" 



37rV2 



-tdf 



(3.10) 



(ii) The limit x ^ 1 

With the progressive loss of energy to the ambient stars, the orbital velocity of 
the secondary BH may increase above a (see §4 for a discussion). For a; S> 1, the 
evolution equations, after considerable manipulations, can be cast into the following 
form 

where in /C(e), and Z{e) we confine the dependence on e (the details are presented 
in Appendix B). Equations (3.11) and (3.12) show that the time evolution of E and 

J is not exponential in this limit but follows approximately a power-law. 

These equations contain the relevant time scales and t> that depend, in this 
limit, on e and on a new quantity 

X., = , (3.13) 

representing the velocity the BH would have as if it were to move along a circular 
orbit; this velocity is expressed in units of the dispersion velocity of the background 
stars. In Xdr we mainly retain the dependence of the time scale on the energy E (see 
eq. [2.18]). 

The resulting time scales are 



/C(e) 



5xl0^(^) ^yrs, (3.14) 
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It is of interest to notice that, in this regime, the time scales of orbital decay and 
rd. are longer than the "reference" time scale tdf (i-e., the scale for x <^ 1), and the 
dilution factor is of ~ x^, as DF on oc (equations [2.6] and [3.1]). Considering 
further their dependence on e through /C and Z, both time scales are found to become 
progressively shorter with increasing eccentricity In Figure 1, we show t> and t> 
versus e and note the decrease of the time scales as e — > 1 : As an example for J, 
T^(e = 0.9) -0.lT^(e = 0). 

The release of energy by the decaying binary is a source of heating for the back- 
ground. For our choice of parameters, however, the relative enhancement of the dis- 
persion velocity is negligible, being Acr/cr ~ G M2/r]rc/N ma"^ ~ M2/MC -C 1. 

3.2 The Criterion for the Growth of the Eccentricity 

In the decay of the orbit, energy and angular momentum losses to the field stars 
may lead either to the growth or the decrease in the eccentricity of the system, the 
effect having important consequences on the life-time of the binary. We here formulate 
two simple and general criterions for the growth of e, the first applies to the motion 
in the harmonic potential and the second to the motion in the Keplerian potential. 

(a) The Harmonic Potential 

From equation (2.13) we can express the eccentricity in the following form: 

TA + rp 



e' 



n J \n 

The time derivative of the eccentricity therefore reads 



de^ 2 ■ 



(3.17) 



It is obvious that the sign of e depends on the competition between the energy loss 
and the angular momentum loss, i.e., on the sign of EJ^EJ. We are therefore led to 
introduce a general criterion for the growth of e using equation (3.17). If the condition 

{EJ-Ej)>Q (3.18) 
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is satisfied along the BH orbital evolution, the eccentricity will rise. 

Substituting into equation (3.18) the expressions of the energy and angular mo- 
mentum losses (eqs. [3.4] and [3.5]) it is straightforward to verify that (e) = 0. The 
eccentricity remains thus frozen to the value set at the outset of the BH evolution. 
This holds if only the lowest order in x is retained in equation (2.5). If we include 
the next leading term in the expansion of g{x) we obtain 

(EJ-Ej) = ^^jn^rl-rlf. (3.19) 



The r.h.s. of equation (3.19) is definite positive. We thus find that the eccentricity 
increases during the first stage of the orbital motion of the secondary BH. 
We can also calculate the rate of increase of e 

{e) = \^^^n\rl-rl) (3.20) 




and the associated time scale 

Th = 5V271TBF 

Note that th exceeds tdf by a factor ~ {rdrp^Y : This indicates that the rise of e is 
not severe in this evolutionary phase as it occurs on a time scale longer than the time 
for energy and angular momentum losses. 

(b) The Keplerian Potential 

In the region rcusp < r < rs, the orbital motion of the secondary BH is described 
by the equations (2.18) and (2.19). They relate e to J and E to give 

J^ = V^(l-e^). (3.22) 
We can thus derive the expression for the time derivative of the eccentricity 

fj 1E\ , , 

Similar to the previous case, we can introduce a general criterion for the growth 
of the e. Using equation p.23| ) the condition is 

0<-^. (3.24) 
(J) J ^ ^ 
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Below, we apply this criterion of growth using the results of § 3.1. 



(i) The limit x < 1 



Equations (3.8) and (3.9) indicate that the eccentricity is constant along the mo- 
tion (i.e., (e) = 0) since {E)/{J) = —2E/J sharply. This result holds if only the 
lowest order in x is retained in equation (2.5). If we include the next leading term in 
the expansion of g{x), the long term changes of the energy and angular momentum 
are found to be characterized by a ratio 



(A 
(j) 



2E 



X 



X 



1 - 



1 - 



3 GM 
207r(72 a 

3 GM 
207r(72 a 



(1 - e2)-V2/, 

(1 - e^y/^h 



(3.25) 



that depends now explicitly on e (the integrals Ii and I2 are given in Appendix B). 
Prom equation (3.18), the criterion for the "growth" of e is fulfilled if 



1 < 



1 72(e) 4 - 3(1 - e2)V2 



(l_e2)i/2 



(3.26) 



This inequality holds for > 0, and equation (3.26) is thus satisfied Ve e [0, 1[. In 
the limit a; <C 1 we thus find that DF on the secondary BH acts so as to increase the 
eccentricity, during orbit decay. The time scale for the growth of e can be estimated 
considering the evolution equation for the eccentricity. 

Following Bertotti and Farinella (1990), we compute the secular change of e (see 
Appendix C). After some calculations, we are able to find a simple analytical expres- 
sion 

8x GM (1-6^)1/2 



l_(l_e2)i/2 



(3.27) 



5(27r)V2cr3 a e 

Note that (e) is positive definite for any value of e, and the characteristic time of 
growth for e is 



tdf 



57rV2 

8 xl, 2(1 - e2)V2 [1 _ (1 _ e2)i/2] 



e 



yrs 



Xc. 



J 2(l-e2)V2[i_(i 



e2)V2] 



yrs. 



(3.28) 



The behavior of r< is reported in Figure 2. t< is found to increase as e — > 1 but the 

deviation is significant only for extreme values of e, very close to unity. This time 
scale is longer than the scale on which energy and angular momentum are lost (eq. 
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[3.10]), and the dilution factor is ~ Hence, only a weak rise of e is expected 

to take place, in the early phase of binary decay. 

(ii) The limit x ^ 1 

In this limit, the ratio {E)/{J) reads 



(J) 



/GM\^/^ /C(e) 
K^J 'Zie) 



2E 



Zie) 



(3.29) 



Considering the dependence of /C(e) and Z[e) (drown in Figure Bl and B2 of Ap- 
pendix B) we find that criterion (3.24) is satisfied Ve. As a result, e increases since 
angular momentum is lost by DF more effectively than energy, as for the opposite 
limit. 

In particular, at small eccentricities (e <S 1) the ratio (£')/( J) takes a simple form 
that we derive in Appendix B: 



(A 
{j) 



-f(l-3e^) 



< - 



2E 
T' 



(3.30) 



From equations (C.l) and (2.6), the long term change of e is governed by the 
following equation (for the details see Appendix C) 



TT 



(3.31) 



where W(e) is a positive function given in Appendix C. 
>V(e) takes a very simple form in the case e -C 1 



W{e) 



1 + Y + '^'"\ 



(3.32) 



The associated time scale is 



W(e) 



4 X 10« f ^ 



W{e) 



yrs. 



(3.33) 
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The behaviour of r> is shown in Figure 3. At small eccentricities the time scale is 

1 _ r ^ ^ _ 1 

(3.34) 



e 1 3 



Contrary to the case for a; ^ 1, the evolution time is found to decrease with increasing 
e : Thus, eccentric orbits become more eccentric on a progressively shorter time scale: 
We may refer to this effect as "runaway" increase of e. 

(iii) X ~ 1 

We have studied so far the long term evolution of E, J and e in the two asymptotic 
limits oi X. To complete our analysis we explore numerically the evolution equations 
(3.1) and (3.2) using the complete expression of g{x). 

The long term evolution of E and J when {x) ~ 1 is found to be intermediate the 
two bracketing limits. Accordingly, we find that the eccentricity increases during orbit 
decay since angular momentum is always lost more effectively than energy. Whereas 
in the hmit (x) 3> 1 the losses at pericenter are weaker than at apocenter (eq.[2.6]) 
the opposite holds when (x) <S 1 (eq.[2.5]). Nonetheless also in the last case, (e) 
maintains the same sign. Pervious arguments based on comparisons between losses 
at different orbital phases were used to derive conclusive statements for the growth 
of e; the above finding thus shows the weakness of that approach. 

We summarize the results of numerical integrations in Figure 4 , where {E) /(J) in 
units of —2E/J is shown as a function of e, for selected values of Xdr- It is remarkable 
that regardless the specific value of x, this ratio maintains below unity over the 
complete domain [0,1 [, and becomes smaller with either increasing eccentricity (close 
to unity) or increasing x. This is in accordance with the behavior of (eqs. [3.28], 
and [3.33]). 



4. POSSIBLE EVOLUTIONARY PATHS FOR THE BH BINARY 



In §3 we have shown that the BH eccentricity increases slowly in the region < r < Vc 
where the gravitational field of the core dominates. As a consequence, about r^, e 
has a value close to the one acquired about Tc. The orbital parameters thus keep 
memory of the previous BH evolution in regions r > r^. The action of DF at radii 
r > Tc where the density of the stellar distribution decreases nearly as r^^ have been 
investigated numerically in PR. In this work it is shown that the BH reaches Vc with 
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a small eccentricity: It follows that the BH trajectory is "prepared" to have e <^ 1 at 

Starting from these considerations, in this section we explore the long term evo- 
lution of the secondary BH and in particular the evolution of the eccentricity e{E, J) 
in the "inner" core region where the "instability" that drives a bound orbit to nearly 
radial infall (e — > 1) may set in. 

Consistency with the model assumptions of a uniform star background requires 
the distance of closest approach not to drop below Tcusp, around Mi (i.e., around the 
focus of the ellipse). In an equivalent way, the apocenter should not exceed tb for 
the motion to be consistent with the Keplerian potential. These constraints lead to 
the following inequalities (for a and e) 

a {1 + e) < T]^/^ Tc (4.1) 



a(l-e)(l + e) 

^ ->rirc. (4.2) 

e 

For sufficiently high values e and binding energy —E, the BH can probe the region 
within the cusp where density gradients are present and where cr ~ v; this region is 
not accessible to our model. 



4.1 Selected Paths 



We begin the numerical integration of the evolution equations (3.1-3.2) at a dis- 
tance tti = 20 pc (for which (xj) ~ 0.3) and choose a value of ?7 ~ 0.05 (following 
BBR). Along the course of the BH evolution a is kept constant (see §3.1). 

The eccentricity is treated as a free parameter; small initial values are preferred 
being consistent with the previous BH history (PR). In the evolution equations, the 
drag by the ambient stars is computed using the complete expression of g{x) (eq.[2.3]). 
The means (eq. [3.3]) are evahiated numerically and a test on the accuracy of the 
calculation is performed using the analytical results obtained in the asymptotic limits. 

Figure 5 shows possible paths for e as a function of the semimajor axis a which 
can be considered in place of the time coordinate. It appears that, independent of the 
initial value of the eccentricity Cj, the growth of e with decreasing a is not severe as 
long as the BH velocity does not exceed the dispersion velocity of the ambient stars 
(i.e, until {x) < 1). At separations where (x) ~ 1, the value of e still depends on 
Cj and the binary keeps memory of the initial condition. The "runaway growth" of 
e occurs only as the mean value of (x) increases just above unity. In Figure 5, the 
function e(a) is tracked also for a < Vcusp (dashed line). In this case the pericenter 
distance rp falls within the cusp radius, and the above inequalities (eqs. [4.1] and 
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[4.2]) are violated. The path below rcusp is shown just to indicate how rapid would 
be the growth of e if the BH were to move in a uniform medium even below Tcusp- 

We now consider an evolutionary path a(t) having = 20 pc, and a small Cj = 
0.03, to follow the evolution of the pericentcr rp(t) = a (1 — e) and apocenter rxit) = 
a (1 + e) distances with time. According to the analysis of §3 we find that in the 
"inner" core region where the BH velocity {v) < a, energy and angular momentum 
are lost to the background stars exponentially, on a time scale tdf ~ lO^yrs. In this 
phase, the semimajor axis decreases exponentially but e does not vary appreciably 
having a longer time scale of growth (oc ty)f/{x)^; cfr. eq. [3.28]): Accordingly 
rp(t) ~ ta(^) ^ (i{t). The binary evolves rapidly through this early phase, until 
(v) ~ a. The dynamical evolution will proceed from this moment onward on a 
different time scale ~ {x)^tj)f ~ lO^yrs, since DF (ex v"'^) weakens as the binary 
hardens. The increase of the stellar density (by a factor ~ 2.5 for a density profile 
~ 7--7/4. ggg Frank & Rees 1976) would however partly compensate this delay in the 
BH evolution. 

In summary, our idealized study shows that we can distinguish two evolutionary 
phases for the BH "binary": (1) an early phase, lasting a few 10^ yrs, in which the 
inspiral of the two BHs occurs along bound orbits of slowly increasing eccentricity 
(the binary is " nearly stable" to the rise of e). (2) a subsequent phase, not described 
by our model, in which the BH probes the cusp region at pericenter, and the uniform 
core at apocenter. The non-local nature of DF prevents us to give a correct estimate of 
the evolution of the eccentricity, during the BH inspiral. According to recent studies, 
the eccentricity would remain nearly constant in the region below the cusp where 
stars can still be treated as a continuum medium (sec PR). Where collisions with 
single stars come into play, e can instead increase again, as indicated in numerical 
experiments by Roos (1981) and Hills (1983). Hence, the evolution of e depends 
sensitively on the detail of the interaction with stars in the innermost region of the 
galactic core, before losses by gravitational waves drive the binary to final plunge. 



5. DYNAMICAL FRICTION AND GRAVITATIONAL WAVES 



In the previous sections we have shown that in a homogeneous isotropic background 
of stars, DF acts so as to bind the binary while progressively rising e. Eventually, 
energy and angular momentum losses by gravitational wave (GW hereafter) emission 
will determine the terminal evolution of the binary. These losses depend sensitively 
on e: For the same energy E, a binary with e ~ 0.9 suffers a loss a factor 1000 larger 
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than for a binary with e ~ (Peters & Mathews 1963). Due to this dependence, 
the transition distance at which the energy loss by DF equals the power emitted 
via GW radiation in a binary is a function of e. at{e) can be estimated setting 
(Edf) — {Egw) where 



GW/ 



32 G^Mf ^. , 



(5.1) 



is the energy loss by GWs (Peters & Mathews 1963); q = M2/M1 is the mass ratio 
and the function /(e) reads 



/(e) = 



1 _L Z3p2 I 37 4 
^ 24^ ^ 96 



'I _ e2)7/2 ■ 

Considering equations (5.1) and (3.11), the equality in the energy losses leads to 



(5.2) 



16 ^ 
5 In A n„m , 



2/11 



2/11 



11 



(5.3) 



Note that at[e) depends on e through /(e) and /C(e) and increases with e. 

The expression of at is consistent with the our treatment of DF only for at{e) > rjrc. 
This inequality translates into a limit on the parameter 77 as a function of e, q and 
of the core parameters: 



V < 



16 G5/2 ^ 
5 In A 



1/2 



m 



.1/2 



(5.4) 



77 depends sensitively on Tc and increases with increasing e. However, for typical 
parameters of a galaxy core, we find that 7] ~ 10^^. In a real scenario DF is effective 
also at smaller radii and the actual value of the transition radius depends sensitively 
on the dissipative mechanisms for energy and angular momentum transfer by the 
ambient stars (see BBR). 

5.1 Braking index 



The simultaneous presence of the dissipative forces of DF and GW reaction close to 
at leads, in this idealized scenario, to the equation for the evolution of the semimajor 
axis a 

64 G3 m V 



5 c> 



-/(e) - a 



21/2^ Em 



(5.5) 



derived using equation (5.1), and (3.11) for x ^ 1. 
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This evolution equation can formally be written as 



4fl + ^a^^/^ 
a 



where 



a = ——M^fxf{e) 

c 



and 



7 



Me). 



(5.6) 

(5.7) 
(5.8) 



7r(GM)3/2 

With use of Kepler's third law, we can relate the evolution of a{t) to the frequency 



f{t) = 1/P of orbital motion of the system: 

Due to orbit decay, we further have from ( |5.6| ) and (|5.9|) 
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a 



+ 7a 



5/2 



We can thus introduce the braking index n 



n 



P 



(5.9) 

(5.10) 
(5.11) 

(5.12) 



that takes the following form 



n 



11 

y 



1 + 2aiV2 



n -1 



a 



11/2' 



,at(e)^ 



-1 



(5.13) 



where now = 11/3 is the value of the index in absence of DF. 

In principle, a measurable deviation from 11/3 would provide an indirect way for 
testing the ambient star medium. In Figure 6 the braking index n (eq. [5.13]) is 
shown as a function of the BH binary separation for various e. It appears that a large 
deviation from 11/3 occurs close to a^, where the two forces are competing. Above at 
the braking index is quite small since DF is the cause of the deceleration of the orbit. 
Below at on the contrary the index is close to new, since GW losses dominate. 
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Recently, Giampieri (1993) considered the deviation from the value new induced 
by mass transfer. The correction term in that case has a sign that depends on the 
mass difference of the two BHs and on the details of mass transfer. The calculation 
of n is here proposed as example for suggesting a way to probe the environment of 
a BH binary. The validity of such approach depends however on the relevance of 
the process considered in the region where an outgoing GW signal is detectable. For 
the case of DF, the GW signal about at would be still very weak and of extremely 
low-frequency. The possibility of revealing GWs about at is thus quite remote. On 
the contrary, mass transfer processes occurring close to the gravitational radius of 
the BHs in the binary, could appear in the measure of n as deviation from the pure 
gravitational effect. 

At low enough BH separations the energy loss to the field stars via single encoun- 
ters lasting less than a period P, might produce a sudden rise in the GW frequency /, 
with a A/// ~ m/M for an impact parameter ~ a. For smaller impact parameters, 
of the order of the gravitational radius Rg{M2) = 2GM2/c^, the " glitch" might have 
a magnitude A/// as large as ~ {m/M)a/ Rg{M2). The potential observability of 
these phenomena would be of considerable importance. 



6. CONCLUSIONS 

We have explored the long term evolution of a black hole pair under the action 
of dynamical friction by a stationary homogeneous and isotropic background of light 
particles. The result applies to the case of a BH pair in the center of a galactic nucleus, 
in regions where the stellar distribution can be approximated as homogeneous. 

Future work should focus on the study of DF in a nonuniform anisotropic star 
background with the aim at exploring the subsequent binary evolution (BBR; PR). 

We gratefully thank B. Bertotti for enlightening discussions and constant encour- 
agement. This research was carried out with the financial support from the Italian 
Ministero dell'Universita e della Ricerca Scientifica e Tecnologica. 
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APPENDIX A 



In this Section we outline the analytical method used to calculate (S) in the case of 
the Keplcrian potential. We proceed along three steps: 

(1) First we specify the instantaneous motion of the secondary BH, described by 
the velocity vector v and the separation vector r. Since DF is a weak perturbing 
force, the motion is determined by the driving acceleration force of the BH primary 
(see § 2). V and r are expressed in term of the orbital phase ■0: They are defined 
uniquely by the instantaneous values of E, J, a and e. For r = r/f we have 

r=-^ T- A-1 

1 + e cos ip 

The velocity is decomposed along the radial and tangential directions to yield the 
following expressions: 



V ■ r — r 



J 



GM 



ail - e2 



1/2 



e sin ijj , 



(A.2) 



Vt 



J_ 



Furthermore, 



1 + e cos ip) 



GM 



GM 



a(l-e2) 



1/2 



^1 + ecos-^) 



(l + 2ecosV' + e^) 



(A.3) 



(A.4) 



a(l -e2) 

(2) Secondly, the perturbing force Vdy is decomposed along three orthonormal 
directions defined by the BHs separation vector r, the velocity v and J. These 
components are indicated with R acting along r, W along J, and T in the orbital 
plane. The three components of the frictional deceleration given by equation (2.1) 
are accordingly to Bertotti & Farinella (1990) (BF hereafter) 



R = 



r 

Vbf • - 

r 



= -X9{x)- 



esinip 



l + 2e cos -0 + 62)1/2 ' 



(A.5) 
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T 



J X r 

Vdf J — 

Jr 

-jX9{x) — 



= -X9{x) 



1 + e cos t/j 



(l + 2ecosV' + e2)V2 ' 



(A.6) 



W 



(A.7) 



In particular, for x -C 1 we have : 



R 



4x 



GM 



3(27r)V2^ 

^5(27r)VV3 
+o{x^) , 



1/2 



a(l-e2) 

GM 
a{l - e2) 



esinip 



3/2 



esin'i/;(l + 2ecos'i/' + e^) 



(A.8) 



T = 



4x 



3V27r(j 

2X 



GM 



nl/2 



a(l - e2 
GM 



5(27r)i/V3 



a(l-e2) 



(1 + ecos'0) 

3/2 

(1 + e cos ■0) (1 + 2e cos + e^) 



(A.9) 



and for x ^ 1: 



R = -2xa 



2«(l-e^) 



esin'^ 



T = -2xa 



GM (l + 2ecosV' + e2)3/2 ' 
2a(l — e^) l + ecosV' 



(A.IO) 



(A.ll) 



GM (l + 2ecosV' + e2)3/2 ■ 
(3) As third step, we express E, J and e in terms of R and T. Following closely 



BF we have 



E — //Vdf • V 

— ^— [i?e sin ■(/' + T(l + e cos '?/')] , 



(A.12) 
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J 



(A.13) 



and 



e = 



P( l-e^ ) 
27ra 



2U/2 



e cos^ ■i/' + 2 cos + e 



(A. 14) 



1 + e cos 

with P the Keplerian period. The derivatives written above are functions of the 
orbital phase iJj. In Appendix B we carry out their mean, over the orbital motion. 
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APPENDIX B 



In this Appendix we sketch the derivation of the equations for the long term evolution 
of E and J. The calculations involve the explicit integration of equations (A. 12) and 

(A. 13) over the orbital phase ip. Below we indicate the main steps. 

Along a Keplerian orbit, the phase ijj is related to the time coordinate by the 
relation 

dt^^il-e^/^- —di:. (B.l) 

27r^ ^ (l + ecos^)2 ^ ^ ^ 

This equation is used to transform the time integral (3.3) to an integral over tp. The 
mean variation of energy E thus reads: 



iix fGM\^/\^ f^^ ,,,(1 + 2e cos V' + 6^)1/2 

{E)^-^i / 9[x{^/J)^———-^-—^d^|J■, B.2 

ZTT \ a / Jo (l + ecos'^j"' 

equivalently for J: 

(j) = _e2)5/2 [^(^)]^ dTjj. (B.3) 

^ ' 27r ^ ^ Jo ^^^^(l + 2ecosV' + e2)V2(i + ecosV')^ 

In the two separate limits x <^ 1 and x ^ 1, the above integrals can be cvahiatcd 
analytically. For arbitrary x{ijj), the calculation is carried out numerically and Figure 
4 summarizes the results of the integration. The analytical expressions in the two 
asymptotic limits were used as a test on numerical accuracy. 

(i) X < 1 

In the calculation that follows, the frictional drag g{x) is given by equation (2.5). 

According to the expression of g{x), we distinguish in {E) the contribution associated 
to the first leading term (oc x) in the expansion of g{x) (indicated with label 1) and 
the subsequent term (oc x^) (indicated with label 3). The frictional drag g{x) is 
known function of ip through x given by equations (2.3) and (A. 4). The resulting 
integral will formally read: 

{E)^{E), + {E)s + o{x^). (B.4) 

As for the energy, the variation of J splits in two terms that will be evaluated sepa- 
rately 

{j) = {j)i + {j)3 + o{x'). (B.5) 
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From equations ( |B.2| ), ( |B.3|) and (2.5) we find: 



(^)3 
(^)l = 



2x GMii 
37r(27r)V2ff ^ 



(l-e^)V2ji(e), 



57r(27r)V2f^3 \ a ) (1 - e2)V2 ' 



1/2 



where 



^ 57r(27r)V2ff3 

2t (l + 2ecos^ + e2l'^ 



/o (1 + e cos'?/')^ 
After some straightforward algebra we have: 

2n 



(l-eVi(e) 
dip , 



/n e 



/i(e) 



hie) = 27T 



- (1 _ e2)3/2 ' 

27r 

^ (1 -e2)i/2 ' 
4 _ 3(1 _ e2)i/2 



Substituting the values of Jq, /i and /a in (^), (CT, ( [BiSl) and (^) 
equation (3.8) for (-E)i, equation (3.9) for (J)i and 

. 2x/i /GM\^ 4 -3(1-6^)^2 

^ 5(27r)V2(x3 1^ a J (1 -e2)V2 ' 



2x/i fG^M'^\ 



1/2 



5(27r)V2cx3 \^ a y 



(1-6^)1/2 



(B.6) 

(B.7) 
(B.8) 

(B.9) 
(B.IO) 

(B.ll) 

(B.12) 

(B.13) 

we obtain 

(B.14) 
(B.15) 



respectively. 



fii) X > 1 



In the limit x ^ 1, equation 

{E) = 



becomes: 

2 / ^ N 1/2 



a 



cm; 



) /C(e) 



(B.16) 
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where 



/C(e) 



ZTT 



We can perform the integration analytically, expanding 
found 

^ ^ L 4 64 ^ 
The behaviour of /C(e) is given in Figure Bl. 

Analogously for (J), from equation (|B.3D we find: 



{j) 



TT GM 



Z{e) 



where 



Z(e) 



[1 - e'y/' 



2n 



-3/2 



Expanding equation ( p.20|) for e <^ 1 we have: 

Z(e) 

The behaviour of Z{e) is given in Figure B2 



l + }le' + —e' + o{e') 
4 64 ^ ' 



(B.17) 



17| ) about e ~ 0; it is 



(B.18) 



(B.19) 
(B.20) 

(B.21) 



APPENDIX C 



Now we derive the formula of the growth of the eccentricity. Substituting in equation 
(A. 14) the equations (A. 5) and (A. 6), and subsequently averaging we find: 



[MM e + cosV" , , 

(1 + 2ecos'?A' + e^)^'^[l + ecosip)^ 



(i) X < 1 



Retaining, as in the case of E and J, terms to order in g[x{ip)], we write (e) in 
the form: 

(e) = (e)i + (e)3 + o(x='). (C.2) 
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Using the asymptotical expression of g{x), we obtain: 

4x 



{eh 



37r(27r)V2^ 
2x fGM 



57ra3(27r)V2 V a 



(1 - e'r'Q,{e) , 
(l-eY'Qiie) 



where 



/■^'^ (e + cos-?/')(l + 2ecosV' + e^)" , , 

Qn{e)= 7— -T^ 

JO [l + ecosipj^ 



After some straightforward algebra we obtain: 

Qo(e) = 

and 



Qi(e) 



l-(l-e2)V2 



(ii) a; > 1 

Substituting equation (2.6) into ( |C.1| ) we obtain: 

(e) = 2""' 

where 



W(e) = -(l-eV^Q-s/2(e 



Expanding W(e) for e ^ 1 we find: 



>V(e) = -e 



1 + -e" + o(e" 

o 

The behaviour of W(e) is given in Figure CI. 



(C.3) 
(C.4) 

(C.5) 
(C.6) 
(C.7) 



(C.8) 
(C.9) 

(C.IO) 
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FIGURE CAPTIONS 



Figure 1. Time scale r> (in units of 10* yrs) describing the orbital decay of the 
secondary BH in the Keplerian potential of the primary, as a function of e for x ^ 1. 
Top curve denotes the time scale for energy loss (eq. [3.14]); bottom curve denotes 
the time scale for angular momentum loss (eq. ([3.15]). 

Figure 2. Time scale for the growth of the eccentricity e, (in units of lO^yrs) as 
a function of e, for a; ^ 1. The secondary BH is in the Keplerian potential of Mi. 

Figure 3. Time scale for the growth of the eccentricity e, (in units of lO^yrs) as 
a function of e, for x » 1. The secondary BH is in the Keplerian potential of Mi. 

Figure 4. {E)/{J) in units of —2E/J as a function of e for selected values of Xdr 
(0.1, 0.5, 1, 10) summarizing the results of the eccentricity increase for the orbital 
evolution of the secondary BH in the Keplerian potential of M^. Dots indicate the 
analytical expression of (E) / {J) as given by eq. [3.25] (top curve) and eq. [3.29] 
(bottom curve). 

Figure 5. The eccentricity evolution for paths with initial semimajor axis = 20 
pc. The curves give e as a function of a (in pc) for initial values of the eccentricity 
ej = 0.05, 0.03, 0.01 starting from the top respectively. The dashed lines indicate 
evolutionary paths with pericenter distance rp < rcusp- 

Figure 6. Braking index function of the semimajor axis a in unit of at{e = 

0), as defined by eq. (5.13) for selected values of e = 0,0.4,0.7 starting from left, 
respectively. 

Figure Bl. /C versus e, as defined by eq. (B.17): dots denote the analytical expansion 
(eq. [B.18]). 

Figure B2. Z versus e, as defined by eq. (B.20): dots denote the analytical expansion 
(eq. [B.21]). 

Figure Cl. W versus e, as defined by eq. (C.9): dots denote the analytical expansion 
(eq. [C.IO]). 
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